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SUMMARY 

Analytical formulation of the integrated potential 
approach for general unsteady supersonic configurations is 
related to numerical solution approaches using an arbitrary 
finite element mesh. Work remains to be done on adequate 
numerical handling of singular integrals, discussed in an 
appendix. 

Limited results on a planar rectangular wing are 
presented. 
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SYMBOLS 


a Speed of sound in free stream 

o 


A 1 ±’ A 12 ’ A 21’ A 22 


C 


fl,f2, f3 
^2’ F 3 " F 4 


F 


G, G s 


J 


2r 


L, L q 

A 


R 

S 

t 

T 

\y 

u 

v 


r]j|i Z,| 

X.Y.Z 


Influence coefficient matrices in Eqn. 47 
Defined by Equation 4la 

Functions defined by Equations 33 through 39 
A typical function 

Defined by Equations 53 and 56 respectively 
Bessel functions of order 2r 
Singular integral 

Defined by Equations 19, 21 and 22. 

A reference length 
Mach Number 

Radius of a characteristic cone 
Hyperbolic radius 
Body Surface 
Time 

Transformation Matrix (Equation 30) 

Velocity components in the ith coordinate direction 
Normal wash 
Flight Speed ' 

Enclosed volume between leading Mach envelope and 
the characteristic surface 

System Coordinates 

Nondimen sional transformed coordinates 
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SYMBOLS (CONT) 


X =(X,Y,Z) 



A position vector 

(M* - 1)* ' 

Characteristic Envelope 
Leading Mach envelope 
Deformation of body surface S 
Defined in Equation 9 

XMo 

Z 3 

Reduced frequency 

Defined by Equations 25 and 26 respectively 
Velocity potential (Equation 2) 

A column vector of velocity potential 
at the corners of an element 
Particular solution of Equation 3 


yj 

CO 

Sh 


Subscripts 


O 


L , u 


D__ ( 
Dt 


JL_( 


dri 
dV 


A 


) 

) 

) 


Defined by Equation 25 
Circular frequency 

Interpolation function (Equation 27 and 28) 
and Notations 

Receiving field point 
Lower and upper limits 

Material derivative of the deformation of the 
body surface 

Normal derivative operator 
Co-normal derivative operator 
Difference between two values. 



INTRODUCTION 


General . - Developments reported in References 1 and 2 
have demonstrated the approach to unsteady supersonic aerodynamic 
determinations based on matrix/consistent element formulations. 
"Consistent aerodynamic element" implies a transition from real 
continuous variables (downwash, displacement, potential, etc.) 
to numerical discrete variables, such that variationally equivalent 
work is done. Additionally, the arbitrary shape and orientation 
of the (triangular) element enables effective matching of 
problem boundaries. These two features were contrasted with 
general Mach box approaches in References 3-7. 

For shuttle-type flight vehicles with comparatively large 
body diameter to wing span ratio, non-linear body interference on 
the lifting surfaces is very significant in determining the 
unsteady forces involved in flutter. In steady flow, thickness 
may be taken into account in an approximate fashion by a "local 
linearization" and replacement of the nonlinear equations by 
locally linear ones with variable coefficients. Mathematical 
justification is discussed in References 8 and 9 . In the case of 
unsteady flow this approach suggests that the fluid motion might be 
calculated with satisfactory accuracy from linearized, unsteady 
equations containing local values of the flow parameters as 
affected by thickness, mean surface gradient, etc. (References 
10 , 11 , 12 .) 

Additionally the downwash/potential integration formulation 
in Reference 2 has disadvantages when considered for use in more 
complex configurations, and other approaches to the integral 
formulation should be considered. 

The present task was to formulate an extension of the 
consistent element approach involving two steps: 

a) a change from a downwash/potential formulation to an 
"integrated potential" formulation (Ref. 13) 

b) the incorporation of arbitrary element attitude 

to encompass more general configurations and space- 
variable flows. 
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This report outlines initial developments for such 
problems including 

(1) choice and evolution of an acceptable analytical 
approach; 

(2) construction of "consistent finite element" 
idealization compatible with structural aspects 
of the problem; 

(3) data handling and computational procedures at a 
pilot program level; 

(4) discussion of future generalization. 

It was assumed that some "best" method of defining the 
steady field exists as a basis for superposing the linearized 
unsteady perturbation (e.g. Ref. 14). 

Mathematical Approach 

Due to the complexity of the linear integral equation 
kernels relating downwash and velocity potential or pressure, 
general analytical solutions for arbitrary Mach numbers, plan- 
forms, and frequency parameters are not available. Consequently, 
several numerical approaches, with various advantages and dis- 
advantages, have been developed. The three most frequently 
employed methods are: 

1. The Collocation Method . - This method uses the inte- 
gral relation between the downwash at a point and the 
loading/acceleration potential assumed on the lifting 
surface lying in its forward-Mach cone region, (e.g. Ref. 
15). The simplicity of this method is that it yields 
directly the loading corresponding to any given 
downwash vector. However, a severe short-coming of 
the method is that discontinuous loading requires 
special choice of the loading functions. Also, for 
high frequency parameters, the loading functions 
become highly oscillatory, requiring very close collo- 
cation points in the chordwise and spanwise directions. 
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The Integrated Downwash Method . - This method yields 
velocity potential at any point as an integral 
relation between downwash and source, dipole or 
quadrupole distributions in the forward Mach cone 
region (Refs. 1 to 7 , 16, 17)- In the case of sub- 
sonic leading edges or in "interference" cases, 
downwash or source strengths in the so-called 
’diaphragm surfaces’ are unknown and must be determined 
from the condition of zero diaphragm pressure differ- 
ence. However, at low supersonic Mach numbers the 
diaphragms may be much larger than the lifting surface 
and the approach becomes very inefficient from compu- 
tational time and storage points of view. Another 
serious drawback of this method lies in difficulties of 
determining the downwash/source strengths along singu- 
lar leading edges. 

For problems with interference the method requires both 
source -downwash and potential-source strength solutions 
adding to the total computational effort needed. 

The Integrated Potential Method . - (Refs. 13, 18, 19) 
This method yields the velocity potential at a point 
as an integral involving the velocity potential taken 
over the forward Mach Cone region. Since the potential 
difference across diaphragms is zero, consideration of 
these regions is not required. 

The advantages of this method are: 

1. diaphragm regions are avoided ; therefore 

problems involving determination of singular 
downwash terms/source strengths in wakes 
are eliminated; 



2. velocity potentials can be more easily 
approximated, since they are continuous; 

3. it should be computationally 

efficient compared with method 2, especially 
for complex configurations. : 

Discretization 

For thick bodies in supersonic flow, local velocity 
and Mach number distribution vary considerably. Any numerical 
method using characteristic grid systems requires change of 
the grid system with local Mach number. In the present approach 
"consistent finite element" formulations are used in a fixed 
grid system which can account approximately for small variations 
of local mean flow parameters; e.g., local Mach number and 
velocity vector due to non-slender configuration. 

The basic advantages of the finite element integrated 
potential formulation are formidable. 

1. Excellent boundary representations of leading, 
trailing and tip edges of wings, etc., compared 
with box and "sub-box" representations; 

2. Possibility of tailoring other elements to 
nacelle and body representations (e.g., ring and 
conical "shell" elements); 

3. Flexibility with respect to element type, size, 
and orientation so that grid point density can be 
varied compatible with spatial loading variations 

4. Kinematically consistent downwash, loading and 
generalized force derivations; 

5. Efficient handling of wake interference in that 
wake edges can in effect, constitute a "single 
element 



6. Feasibility of using common structural-aero- 
dynamic grids is convenient for optimization 
problems ; 

7. Formulation for locally varying Mach number 
is possible, given a steady state flow field 
distribution. 

These advantages are expected to be even stronger 
in more complex geometries such as wing-body interaction in 
terms of feasibility of adequate modelling and economy of compu- 
tational effort. 


INTEGRAL EQUATION FORMULATION 


The following development parallels that of Ref. 13 to 
a large extent, the basic differences being (i) the extension of 
the formulation from planar surfaces to arbitrary bodies (within 
the bounds of a linearizable approach), and (ii) the retention 
of the local Mach number parameter y9 so that some approximate 
handling of its effect in a variable Mach number field is 
facilitated. 


The 

satisfies 


velocity potential $ (such that M. = 

U i ) 

dXi 

aV . = a? 

c^t 2 <5cc,dfc doc., 2 


which transforms , using 

0 = (p eaop (- iX X ) eogofiATl 

P 

tfcp L(9)» o 

iY 2 bZ z dx z Y T 


(1) 

( 2 ) 

(3) 


in which x _ -x-'/pir , y = u,/i r , z _ Vc are the 

non-dimensional coordinates with respect to the reference length £}• 
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Defining X 0 = (X 0 , Y 0 , £„) a field point, and 

X. =(X 0 _X), a = (Y<,-Y), Z = (z 0 -r), 

the characteristic cone from 


X D is the surface 


x 2 - y 2 


, z — 


= x a — r 2 - = O 


(M 


If y is some particular solution of (3), the generalized 
Green's theorem states 


J [> L (<p) _ (p L (>)] d V = J _ v |£] d S 


(5) 


In equation (5) & i* _ l _ m — 

3v ~ ax 3Y 


4* 


n 


b Z 


is a differential operator, S> is the co-normal such that 

t'-Yn.n. are direction cosines of the normal on S into V. (Fig. 1) 
S and V are taken as the surface and volume defined by 
the leading wave front, the forward characteristic surface 

from and the body surface. Since L (<p) = L (y) - 0, 

Equation (5) becomes 

_ > <L£1 d S — O 
3 v J 

a 

On the wavefront P u , cp and = O , 

is chosen such that yj and P are zero on P, 

^ 3 v 

becomes 

f ( <f _ y/ 3 f\ d S _ O 

J e.j 3 ° v 

Ref. 13 gives the appropriate Y*'* as 


si* 


( 6 ) 

hence if V 
Eqn (6) 

(7) 


r t = Sip ^ e ( 0 4 1 ) 




( 8 ) 
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( 9 ) 


(which are singular on ac — t ) , and where 

q? = k 2 ( jc z _ r 2 ) , e = * 

v ' r- 

By the common potential theory maneuver of excluding the 
singularities by . small limiting surfaces and proceeding to a 
limit, Eq. (7) becomes 




/cR 


where 


2-rrJ l R 


\ f co sUR \-ZV]AS 

) \ r / as> . 


R 2 = oc 2 _t 2 

Using Eq. (2), and dropping exp ( LX"T ) throughout 

( coskR ) + LkeM 0 («o|kR)^(X) 

^cos k < 5 $ 1 k M.. f X — X))dS 

Defining the *normal" derivative at a field point X a 


( 10 ) 

( 11 ) 


(12) 




+ m 'o 4- + *» 


"^« 0 u ax 

the normal velocity at this point is 




n, = ii 

° e r a 7t 


Then, 


a0(x o ) 


^7. 


ur — _ j_ a£(x„) 

^ <r 

r J (o 1 


+ 


l k M 0 1 coskR ] d> (X) _ cos UR (X) clS 

R J R a-*> J 


(13) 


(14) 


(15) 
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2£(Xo) = 


= ibJ^ ( - ikMo(x °- x) [{^( t2 i kE ) + 

k 2 Mj l, 0 L cesk R + Lk M 0 p| ( cosk.R ), < 

_ (<£lb.%) _ t k M 0 £, o cos k R j ^ (X)j ^ S ( 16 ) 


For panels on thin wings and on bodies (i.e., with flow 
on both sides and one side respectively) handling of the A0 
(difference) and aM requires some comment, condensed in the 
following: 




Panel of a 
thin surface 

Nonzero if 

interference 

occurs 

Nonzero in interfering 
case 

Zero for non-interfering 
case . 

Panel of a 
body 

Nonzero 

Nonzero in general 


Allen & Sadler, (Ref. 20), Woodcock (Refs. 18, 21) con- 
sider d */av to be identically zero for plate-like elements. 
However, is not zero for interfering cones such as 

will occur in wing-body interference. This term can and will be 
considered in future developments in this work. 


COMPUTATIONAL METHOD 
General Formulation 

The integral equations (12&16 ) for the potential are so 
complex that closed form solutions are not considered. Numerical 
integration methods are invariably employed for which Refs. 20 & 

21 for example, use characteristic grid systems. In these analyses. 


12 








non-singular terms within each element are expressed by means 
of interpolation functions while singular terms are evaluated 
considering only the finite part of the integrals. For planar 
elements, the singular terms are few in number and the singular 
integrations may not affect the results significantly. How- 
ever, in non-planar cases the integrals in Equation (16) are 
strongly singular and there are more terms than in planar cases. 
This section discusses the possible means of reducing the 
number of singular terms, and the extent or domain of singular 
functions . 


The kernel ( cosk in Equation (16) exhibits square root 
singularity as the characteristic surface is approached. This 
form of singularity can be expressed by a well-behaved Bessel 
series Refs. (3) and (22). Performing the required normal deri- 
vation and simplifying the terms, the velocity potential and 
downwash integrals can be written as 


J- 


1 


eocp[-^kM 0 &xj|j . |^(|^) 

_ ailJS 

ay 


av 


]* 


(17) 


and 


the downwash integral as 


1 

£irt r 



k%MX£- 


dS (!8) 
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In these 

equations. 


Mx) 9 

.(_LL„ + 3L.(X)L Qi) + LkM„|hLCX)_6u(X)jj 

(19) 

in which 

LLo = Ckot 4 4 

(20) 


L(X)= -l AX + m AY + n AZ 

( 21 ) 


L o 0Q= L AX + n, 0 AY + n 0 A:Z. 

(22) 

and 

AX =(X„_X), AY= (Y 0 _Y), AZ = (Z„_Z) 



Watson's relation- (Ref. 22), is 

cos k R ^ yj 

R ~ “ a Y (23) 

in which 

4^)= X(fv*)sin ‘ ’(I) + Sin(2f5io _1 (i)) (24) 

with X = Yo-Y ( 25 ) 

and yU. = >J(Xo — X) — (Z. o ( 26 ) 


The velocity potential in Equation ' (IT) is singular only 
at the apex of the characteristic cone, while the normal wash 
integral has singular terms at the apex and on the characteristic 
surface for a general configuration. However, for a planar case 
the singularity in the normal wash term is limited to the apex 
only. 

The numerical integrations of Equations (17) and (18) can 
be performed on any convenient discrete element basis such as 
triangular, quadrilateral or conical ring elements. This ideali- 
zation requires description of the doublet (A 9 ?) or source 
distributions within each element in terms of its nodal or boun- 
dary values. These distributions might be limited to linear 
or quadratic expressions in space coordinates. 
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Triangular Element. - These are the simplest elements 
suitable for many irregular configurations. Curved surfaces can 
be represented adequately with a sufficient number of flat 
triangles . 

Quadrilateral Element. - Consideration of quadrilateral 
elements in place of triangular elements is valuable in terms 
of computational efficiency for a given flat lifting surface. 

Conical Ring Element.. - For axisymmetric bodies of revo- 
lution it is possible to use a circumferential description of 
normal wash and the corresponding doublet/source distributions 
in the form of Fourier sine or cosine series. 

To demonstrate the method consider the triangular element 
with a linear distribution of doublet or source strength, i.e.. 




and A0 and A are the values of the source and doublet 

strengths respectively at the vertices of the triangle. Using 
( 27 ) and (28) in Equations (17) and (18), the velocity potential 
and normal wash integrals for each element in the domain of 
influence are given by: 
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and 


TjlT — 


_ c 


Znl r 


expf_ 

% 


ikM„AXj 


(F 3 feyA Fa avl a 4 

Lax axVsr/ + aYJ 


_(n_ a_ /av>\ , r; a si Ids 
(ax ^x^y/ ay J a=? 


(32) 


respectively. 


The following definitions are observed in equations (31) 


and (32); 

i « Uxyn.(x,Y) (33) 

j , M<M.m(X,Y) (34) 

f a = Yl(X.Y) (35 ) 

F; = y(X)Yl(X,Y) (36) 

f 2 = k\^u - MX)L(X) ) n, (X, Y) ( 37 ) 

R 

F 3 = L„CX)n,(X.Y) ( 38 ) 

£ = (.kM^flCX.Y) (39) 


The general forms of the integrals in. Equations (31) and 
(32) are of only two types, i . e . , involving single and cross- 
derivatives of the function V • The function y is well behaved 
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in the fore Mach Cone region, but its derivatives are non- 
analytic on the characteristic surface. The functions such as 
f^ to f^ and F 1 to F^ (with the exception of the term 

l.ooloo/r*- in F^ and F^ which will be discussed later 
under singular integrals) are also analytic throughout the domain 
of integration. 


To take advantage of the well-conditioned function 
it is necessary to cast the integrands in a particular form for 
numerical integration. Consider for example an integral of the 
form 

— L.f^.)ds fU 

ax axvay/ (40) 

in which the function F is analytic throughout the domain of 
dependence. 


exp|_ ik AXj 


In performing numerical integration the infinitesimal 
surface area as can be treated as that projected onto a plane 
whose unit normal is parallel to the normal of the element 
under consideration, e.g., if the normal of the element is close 
to the Z-axis then the infinitismal area 


ds = fW^W dXdY 

= C(x,r) dXdY ( 4la ) 

On the other hand, if the unit normal is close to the Y axis, then 

ds - 1' « 


(4lb) 


Using (4la) the integral is 

.X u , _ r Yu(X) 


I = 
XJ 


exp f— ik M 0 AX 1 
AX 


f dYdX 

ax W/ 


YM 


(42) 


The integrand under the Y integral can be rewritten as 


f a_ ( avA _ 
ax 'ay J ~ 


a_ 

ax 


(sL-C r *Y'i) ~ /IT vA - d far yA 
\ay J ax Uy / aYVaX ) 


+ iiL- V 

+ ax ay 


(**3) 
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Using Equation (43) in Equation (42), performing integration by parts, and using 
Leibnitz's rule for differentiation under the integral sign, the integral become 
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* , r X (X) Y(X) 

J_ \ ML.-^dr , ax ydY 

AX 2 ] 3Y + 5>T§Y'^ 

L X(X) \(x) 



The first and the last terms' vanish for elements with no 
sides parallel to the Y-axis. The numerical integration 
utilizes Gaussian weights and pivotal points. Limits of 
integration for complete and partial elements are shown in 
Figures 2 and 3- 


The second set of integrals is the one with a single 
derivative of the function *)Cr i.e.. 

x„ -Yu(x) 

F(X,Y)SV: 3YdX 

aY 

Yl(X) 

Performing partial integration. Equation (44) can be written 

as X, 


If- 


ex.pf— tk.M 0 AXj 


I s 


= f “f'}r(x,r u (X))KX,Y 0 cx))-F(x,x(x))>(x,Y l Cx))l 

J x f XW 

f NET . ... 


jc/drl dx 

£Y 

Y.(X) 


U5) 


(46) 


Using the form of the integrals given by equations (44) and (46), 
the normal-wash integration in Equation (32) will be evaluated 
for all elements in the fore Mach Cone zone. 

Considering all nodes on the lifting surface the velocity 
potential and normal wash in matrix notation are given by: 


and 


A$ = A„ . A IZ A<$ 

I>? _ A*,A<L§- . A 22 l A<£ 






(^7) 
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Woodcock and York (Ref. 21) consider A =>o 

for antisymmetric distribution of velocity potential on the 
upper and lower surfaces of a panel. However, in the presence 
of the interference, since ¥ 0 _ , it appears that a differ- 

ential source distribution is required to satisfy the normal 
wash condition in the second set of equations in ( 47 )- These 
two sets of distributions can be solved from the two sets of 
equations in ( 47 ). However, based on other analyses which omit 
these terms, (Ref. 21 ) the neglect of the term is 

obviously sometimes justifiable. An independent study might 
be conducted to reveal its general significance. 

Influence of the Wake 


Determination of velocity potential for wings with sub- 
sonic trailing edges or tandem surfaces is very much influenced 
by velocity components in the wake region of the forward wing. 

It is a difficult task to determine the equilibrium position of 
a wake sheet even in the case of steady flow for a simple con- 
figuration. However, linearized computation methods assume that 
the wake sheet remains planar with the trailing edge (planes) 
panels. Having defined the assumed position of the wake sheet 
the velocity potential difference across the sheet can be evalu- 
ated in terms of the trailing edge values only. Since the 
pressure difference across the wake sheet is zero the potential 
difference is given by 


A <l> (X,Y) = A (Y(X» «c P L k Mo ($.)(X_X TE )} 


( 48 ) 


The potential in the wake is therefore dependent only 
on the trailing edge nodal values and there is no need to distri- 
bute elements in the wake region, (in contrast to the integrated 
downwash approach) . Instead line elements of trailing edge can 


be used to describe the variation of A# 


TE 


The integrals involved 


in the normal wash are essentially the same as in Equation (18), 
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with the following terms replaced, i.e.. 


€xp j-t. k M 0 AX j replaced by exp 1 - i kM 0 (X 0 )X 

^ 1 SL „ M 0 

- £ *te] 

Mo J 

and 

iQi _ f X Y llT replaced by tQ* — Y — .Y.. 1 

. ur [y u -y y u -yj 


(^ 9 ) 

(50) 


The upper limits of X are determined from the vanishing 
of the hyperbolic radii of the characteristic envelope. (See 
Figure 4 ) . 


Singular Integrals 

Equation (32) contains singular integrands as the 
integration limits approach the characteristic cone. Accuracy 
of solutions depends critically on numerical treatment of these 
singular functions. For this reason, the number of singular 
integrals involved has been reduced by certain transformations. 
Numerical computation of singular integrals is discussed in this 
section. 


The general form of singular integrals is given by 

-Xu Y (X) 

1 5 - f expf-ik Mo AX] rCX.Y) ^_f^\dYdX 

J ax ax \3Y / 

x <- YJ.X) (51) 


where 


Let 


f(x,y) = («£(X) + L,(x)Lon j 

-Y„(X0 

G(X0 = j ^ (X) + l.QQl(x) 1 !_.(!>') dY 
Xcxo r< 


(52) 


(53) 
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Figure 4 Wake Element In The Integrated Potential Method 
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The first term in Equation (53) is non-singular and 
hence can be integrated by parts. The second term is zero for 
coplanar elements, while for non-coplanar elements it becomes 
singular as the field point approaches the characteristic 
surface . 

Figure (3) shows two forms of element inter- 
ception by the Mach Cone, i.e., 

(i) when the element is away from the receiving point; 

(ii) when the receiving point lies on a node of an 
element . 


Case i 


Put 


G s (Xi) = 


XCxiV 




U(X)L(X) / 3vA dY 

R* SxVSyV 


(54) 


where X^ is the ith Gaussian integration pivot in the X direction. 


Since the function is well behaved over the area A BCD, Fig. 3, it can 
be expressed in terms of its nodal or boundary values, i.e. 


QyX L Y J 

• (55) 

Then performing the required differentiation. Equation (5^) can be 
written as v , . 


G s (Xi) = 


F(XiY) dY 

R* 


* ' Y t (xO (56) 

Let Y ) be the last Y coordinate in the Gaussian integration. 

Then Hadamard ' s - finite part of the integral is given by 

.Yo(Xi) 


G s( X 0 = 


.Yu (Mi 


Y l (X 0 


(F(X it Y)-F(x l) Y*)hY ^ r(X;,Y*) 

R2 


y l cxo 


^ (57) 
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The last .term in Equation (57) is given by 

f U dv - L h/(Ax*-Az*)' 4 , ay \i Y “ Y u 


J r 'uv'h; f a / 

dY = __J U 

R z £[Axf-Az?J'*l V 

Y l (Xj) . 


■w 


Uax^-az 1 )' 4 - ay/ 

(-for AX 2 -AZ 2 f tO) 


\ “i * »u 

'•* Y=Y t 


L JYL (58) 

Considering only the finite part of Equation 58 , the finite part 
of the singular integral Equation ( 57 ) is given by 

G s (x ; ) = f " F(Xj.Y)- F(X;.Y*m Y _ FjXi.Y) ^ r[AX»-AZ , J*- AY' 
J Y ,„. R a £[AX s -AZ.-J ,& [[AX'-AZ-h + ay 


0 - ' 


Y.(*0 


(j-or AX z -AZ 2 o) 


Y 0 (Xi) 


dY ^ FCXi.Y*) 




[ J— — 1—1 

Lay u ayj 


(59) 


Case ii 


■ In this case, a singularity exists in the X-direction 
as the integration point approaches . the apex of the characteristic 
cone. Rewriting the. singular, integral using the definition of 
the inner integral given by Equation (5-4); 


r* 4 

u = 

Jv 


exp < ~ t . k M a AX| G(X) <J X 
AX 


( 60 ) 


Let X* be the last X'Coordina te in the Gaussian integration then 
the finite part of this integral is 

X u 




£expfi kM 0 AX} 6(X)-. e xp{-ikMoAX*jG(X*)] dX 


+ eXpC-UMoAX' 6 ] G(X*)& AX t 


( 61 ) 


Thus, all integrals involved have been formulated ..for numerical 
computation with a minimum of singular behavior. 
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DISCUSSION 


The first objective of the study, the extension of the 
integrated potential method of analysis to the three-dimensional 
flow-field problem formulation, has been illustrated. In this 
formulation, the varying Mach number-velocity vector is not 
treated explicitly. Rather, such variations will be considered 
in a piecewise fashion, from element to element, and in the 
determination in a sequential manner of the effective domains 
of influence and dependence of nodal source and field locations. 

Some developments in the formulation of the integral 
equation are discussed resulting in improved handling of singu- 
larities on the Mach cone. These aspects of the formulation are 
however incomplete (See Appendix). 

In particular, the integral equations given in Equation 
18 are more complex than those encountered in the planar case. 
They were therefore checked in a simple configuration in prefer- 
ence to tackling the arbitrarily oriented element. Indeed, 
general element codings were outside the scope of the task and 
limited aims have been achieved. In particular, a quadrilateral/ 
triangular element was formulated and coded for restricted circum- 
stances only, with the original intention of. applying it to a 
simple wing-body interference problem. 

Instead, a rectangular wing of aspect ratio 2 was treated 
and is compared in Table 1 with Reference 23. Agreement is 
limited and improvement is desirable. 

After some investigation the source of inaccuracy (suspect- 
ed as program coding for some time ) was found to be in the formu- 
lation and computational treatment of the improper integrals. 

The present status is discussed in the Appendix. 
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TABLE 1. 

COMPARISON OF GENERALIZED AERODYNAMIC FORCES FOR RECTANGULAR 

WING OF ASPECT RATIO 2.0 AT M =1.2 

o 

Mode 1 - Heave 

Mode 2 - Pitch about mid-chord 


Method 

Present 

Reference 23 

Number of 

• 




Elements on 


5 


30 


Chord 





k 

G.F. 

Re 

Im 

Re 

Im 

0 

Q 11 

3-31 


3-75 



Q 22 

-0.5 


-O.38 



Q 11 

0.0644 

0.938 

0. 188 

1.009 


Q 12 

3.176 

0.027 

3.369 

-0.499 

o.3 

Q 21 

-0.043 

-0.147 

-0.0045 

-0.130 


q 22 

-0.4685 

0.273 

-0.4114 

0.164 


Q 11 

0. 068 

1.730 

0.348 

1.640 


.<>12 

3.057 

0.3034 

2.838 

-0.294 

o.6 

Q 21 

-0. 166 

-0.107 

-0.2424 

-0.2860 


Q 22 

-0.3307 

-0.503 

-0.379 

0.450 




APPENDIX 


This section discusses alleviation or elimination of 
some of the strong singular integrals encountered in the 
downwash functions (Equation 15)- This requires proper 
ordering of integrations and differentiations. 


Consider the downwash integral given by Equation (15) 


x,Y„ 

^ f eacpp* ik M c AX} — ^coskR ^ lkM 0 P c °s k R | <^(X) 

" r X.Y t 

_ |!^j dXdY 


(Al) 

For flat elements (triangular or quadrilateral) the 
normal and co-normal derivatives can be expressed in terms of 
local normal and co-normal derivatives, i.e.. 


and 


a_ = . 

ay ax 

a e a 

= 4- t-oi i— 4- 

ax 0 



(A2) 

(A3) 


using Y and Z component velocities expressed in terms of the 
local Z 7 component (See Figure Al). 

Using the Bessel series expansion for ( CoskR)/R 
(See Equation 24) the expression in Equation (Al) may be 
written as 




W _ 1 


dn' 2n£ r ar> 0 ' J 


expi-ckH,AX}.Jp|_(^ (M) 

f 

+r>'i- PvV, l kM 0 e^v|<6(x)_^v wW' 

dz'\ay'; + ay'J r ay'av'J 
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'O' 


, -T are the system 

coordinates 

r ,1-,' are element loc 

X > Y » 2 ' ‘ coordinates g 

2 * normal to 
plane of the 

element 




. „ and Componer 

Local Coordinate Eleroen t 

Cosines of a UP 


n ta of Direction 


APPENDIX (CONT) 


The last two terms can be Integrated with respect to Y 

/ 

in terms of yj at the upper and lower limits of Y(X), 

Should the limits be on the Mach line, the Sin -1 (%) 


term become 


s ± 7r/z and Y'=s>±!' JoCk/*) , the derivatives of 
which exist. 

The singular nature of the downwash integral stems from 
the first two terms i.e., the integrals of the form 

I, - ^ f feacp i-okMoAX} 0(K,Y) dXdY' 

JJ dX VdY / 


(A5) 


and 


.. |_ J|eip|-il<M 0 AX} j6(X,Y) dXdY' 


(A6) 


Handling of these two integrals will be typified in this section 
by simple examples. For the sake of simplicity let £ = 0, 

(representing a planar wing case). Then the downwash integral is 
represented by 


-ur _ 1 


dZ' j 


X„Y„'(x) 

exphikMoAX} S (|V) 5&(*,Y')dXdY' 

X,Y l (X) 


Since Y>) and Yu(x) are independent of zJ } integration 
with respect to Y' results in 


(A 7) 


U - 


r y | Yu N r s (o 

J | _ i £_J y 6Y J j dx 


2rxl r d Zlj v' 

^° J X U Y«. ’ ‘'YJ 

(since, in a triangular element,^ is linear in X,Yf) 


(AS) 
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APPENDIX (CONT). 


At the upper and/or lower limits of Y^on the Mach line, the 
Y> and Jy/ terms are purely functions of J 0 (k>«.) > the 
derivatives of which exist. Hence the normal derivatives 

^ or jL. can be performed under the integral sign. 

8n<; 2Z' a 

/ 

For illustrative purpose let 0 Cx»Y) be constant 
within a triangle. Then Equation (A8) may be written as 


X 

•ur_ <t> £ f ea.phikM.AX} fyfrxVKXXY] dX 

“ 2nt r azj Sz' L J 

■Xu 


(A9) 


where 


^(xX) = V=+/V ar 
= ±|j 0 (i<^) 


for the sides of a triangle not 
on Mach line. 

for a side on the Mach line, in which 


sin'll) 

* r = (2 r sin- '(^ 

* r 


For a diamond element formed by two triangles (Fig. A2), the 
contribution from lines CB and C'B is 


Xe 

bX. — 0a Idt, . f exp k Mo AX | ^ — . ^1 — ( Jo CM)dX 

*c 

= _kf & dx 

Jy k/A* 

A c 

where y-c * = Xo "X 

Note the regularity of the function, kt. ld^A - i 

yU-^O AyU. 


(A 10) 
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which normally would have been singular. For the steady 
state case (k«o) the contribution to the downwash with 
constant (p in element 2 is zero. 

► Y 

/ 

/ 

/ 

/- 


Figure A2 

♦ ~ ' 

The contribution from the triangle (1) can be obtained by 
differentiation of the and functions with respect to 

Z and Z e and taking the limits as Z.. Since this 

is a non-singular element, no elaboration is necessary. 

A similar treatment can be applied to elements inclined 
to the flow direction (i.e., with L^O ). 

This illustrates how the strong singular integrals 
appearing in Equations (12) and (15) can be reformulated so 
that errors are minimized in the numerical stage. 
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